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Abstract 

In the first part of this contribution we prove the global existence and uniqueness 
of a trajectory that globally converges to the minimizer of the Gross-Pitaevskii en- 
ergy functional for a large class of external potentials. Using the method of Sobolev 
gradients we can provide an explicit construction of this minimizing sequence. 

In the second part we numerically apply these results to a specific realization of 
the external potential and illustrate the main benefits of the method of Sobolev 
gradients, which are high numerical stability and rapid convergence towards the 
minimizer. 
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Introduction 



Ever since Bose-Einstein condensation was realized in dilute bosonic gases in 
1995 [TI2] . experimentalists as well as theoreticians have been interested in de- 
scribing the experimental results with an accurate yet also efficient theoretical 
framework. The most common approach in this context is the Gross-Pitaevskii 
equation pilfS] which dates back to the 1960s and provides a meaningful de- 
scription in the regime of weakly interacting bosons. Recently, it was also 
shown rigorously that the exact ground state of dilute, interacting bosons 
in an external trap can exactly be described by the minimizer of the Gross- 
Pitaevskii energy functional in the appropriate asymptotic limit [6]. 

Even today it is of utmost importance to have an efficient way of finding the 
ground state of weakly interacting dilute bosons in an arbitrary external trap 
as this state is either of interest in itself or the starting point for a time evolu- 
tion. There exists a large number of numerical approaches to find this ground 
state, among them the density-matrix renormalization group method [7], the 
Multi-Configuration Time Dependent Hartree method [5|^|TU] . the Numerov 
method [llj, the fast Fourier transform method p!2] using steepest descent 
or imaginary time propagation, the finite element method or the method 
of Sobolev gradients [H]. However they differ significantly in their numerical 
stability and efficiency and the extension to three-dimensional problems is of- 
ten accompanied by technical difficulties or problems with the efficiency. In 
this contribution we present the method of Sobolev gradients in detail, since it 
provides a very powerful framework for analytical investigations and addition- 
ally offers the possibility of combining analytical as well as numerical results 
in a natural way. This method has before been applied, e.g. in the works of 
P!^16|17j . yet our study differs from these previous works, since we obtain 
global convergence for the trajectory that leads to the minimizer. Using our 
finite differencing scheme, we are thus guaranteed convergence as the mesh 
goes to zero. Furthermore our results are applicable to a very general group 
of trapping potentials, both numerically and theoretically. 

We first discuss the relevant properties of the Gross-Pitaevskii energy func- 
tional which allow us to prove the global convergence of the minimizer of this 
functional with the help of Sobolev gradients. Having shown this result, our 
numerical simulations are guaranteed to provide the correct ground state in 
the limit of an infinitely fine grid on which the Gross-Pitaevskii equation is 
solved. Additionally, we highlight the numerical advantages of the method 
of Sobolev gradients in the second part of this contribution, where we discuss 
one-, two- and three-dimensional simulations to obtain the ground state of the 
Gross-Pitaevskii energy functional in the presence of a generalized Mexican 
hat potential. 
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To set the stage we begin with the definition of the (i-dimensional Gross- 
Pitaevskii energy functional E{iIj) for N interacting bosons of mass m in an 
external trapping potential Vtrap{x) 




where g denotes the coupling constant which for d = 3 reads g = Anti^as/m 
and is determined by the s-wave scattering length a^. The wave function tlj{x) 
is complex valued on R'^ and we are interested in minimizing the energy func- 
tional subject to the normalization condition 



We restrict ourselves to external potentials that are measurable, locally bound- 
ed and tend to infinity for \x\ —>■ oo. Thus the potential is bounded from below 
and without loss of generality we consider the minimum to be zero. Neglecting 
the quartic self-interaction term in the energy functional, the minimizer is 
equivalent to the ground state of the standard Schrodinger Hamiltonian 



H^- — A + Vtrap{x) 

2m 

which provides a natural energy unit of hw. For the most commonly used 
example of an isotropic harmonic oscillator Vfrapix) — miu'^x^/2 this can be 
shown trivially and the corresponding ground state wave function is a Gaus- 
sian. From the energy unit we can derive an appropriate length unit which is 



given by = \Jh/muj. By rescaling the original energy functional in terms of 
these units, that means x = GqX and E{ip) = hiuE^u), we obtain 






IRd V2 





where u{x) = ao'^^'^ip{x), Vtrap{x) = Vtrapix) / fku and g = ao '^g. The normal- 
ization condition remains unchanged and reads 




(2) 
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For the following discussion of the minimization of this functional, we restrict 
ourselves to a bounded domain D G M!^ and minimize the energy functional 
([T]) subject to the constraint ([2]). This problem has already been investigated 
in many studies, both theoretical and numerical ( |18f 19ll2Uf21f22ll23f24] and 
references therein). The existing result for convergence is the existence of a 
minimizing sequence that converges strongly in the norm to a member 
of the Sobolev space H = H^''^{D,(C). In contrast, we will show via steep- 
est descent with a Sobolev gradient the strong convergence in the H norm. 
This is our main theoretical result and the convergence is obtained for a wide 
class of trapping potentials. Furthermore, we do not only prove the conver- 
gence in an abstract sense but actually derive the minimizing sequence z(t) 
in a constructive way, which allows us to use this sequence for our numerical 
simulations. 

The method of using Sobolev gradients to obtain stationary solutions of an 
energy functional provides various advantages compared to other methods 
which are currently being used. First, since the convergence is in the H norm 
we not only get that z{t) converges in the norm to a member of H, but we 
also get that the partial derivatives di{z{t)) converge in the norm, for i = 
1,2, ... ,d. Not only is this a stronger form of convergence than what has been 
proved so far, but it also demonstrates that the numerical simulations based 
on this method will be well behaved in two ways. Since the continuous version 
of steepest descent converges, it is clear that for an increasing number of 
discretization points in our numerical simulations, the solution of the discrete 
minimization problem converges to a member of H. Additionally, one can 
expect that the convergence of the discrete problem will be smooth since the 
convergence of the continuous problem is in H. We outline these benefits 
of the Sobolev gradient in more detail in the numerical part. Eventually it 
is the combination of our theoretical results and their practical numerical 
implementation which constitutes the main new contribution of this paper to 
the problem of finding and analyzing critical points of the Gross-Pitaevskii 
energy functional. 



1 Analytical results 

The goal of the analytical part of our contribution is the proof of the global ex- 
istence and uniqueness of a trajectory that globally converges to the minimizer 
of the Gross-Pitaevskii energy functional. To arrive at this goal, section 11.11 
provides basic properties of this functional which are essential for our proof. In 
section [L2] we define the gradient that we will use for the minimization and de- 
scribe how we incorporate the normalization constraint of the desired ground 
state into this gradient using a particular class of projections. We set up a 
trajectory using this gradient and consider the resulting evolution equation. 
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Of the final two sections of the analytical part, section [L3] deals with proper- 
ties of the class of projections we obtained. We will need these properties in 
section [13] in order to obtain global existence and uniqueness as well as global 
convergence for the minimizing trajectory. 

To avoid confusion, we define the notion of a Frechet derivative before we 
start with the analytical details, since it will commonly be used throughout 
this paper. If G is a function from one Hilbert space to another, we denote the 
first Frechet derivative of G at m by G'{u) and by G"{u) we denote the second 
Frechet derivative of G at m, provided these derivatives exist. Furthermore, we 
remind the reader of the definition of Frechet differentiability. 

Definition 1. Let G : Hi ^ H2 be a map between the two Hilbert spaces 
Hi, H2. G is Frechet differentiable at u E Hi if there is a continuous linear 
transformation : Hi H2 so that for all e > there exists 6 > so that 
for all h E Hi with \\h\\Hi < S the following inequality holds 

\E{u + h) -E{u) ~Tuh\ 
II^IIhi 

We write = G'{u) and G is twice Frechet differentiable if the map u — >■ 
E'{u)h from Hi to H2 is Frechet differentiable. 

In cases where G is known to be Frechet differentiable, the Frechet derivative 
and the Gateaux derivative of G coincide. The Frechet derivative can thus be 
computed in the following way 

GT.)ft^lim^<" + "'>-^'"' 

' t^o t 

and is similar to a functional derivative which is commonly used in Quantum 
Field Theory. 

1.1 Basic properties of the energy functional 

In this section we discuss the basic properties of the functional ([1]) which are 
necessary to obtain our analytical results, in particular the size and shape of 
the domain D of the energy functional as well as its convexity. Furthermore, 
we state the assumptions we make on the trapping potential and show that 
they do not lead to a loss of generality. 

We note that the functional E{u) is well defined for potentials that are non- 
negative and bounded on D, which follows from the Sobolev embedding the- 
orem. Here we assume that D is a bounded domain in R'^ with a regular 
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boundary, that is it satisfies the cone condition. In practice, D is usually de- 
termined by the condition that the ground state effectively vanishes at the 
border of the domain. In the case of an isotropic harmonic trap as external 
potential, that is Vtrap{x) = we can identify the domain to be a ball of 

radius R. This radius is a multiple of the width of the Gaussian in the case 
of weakly interacting bosons, whereas it is a multiple of the Thomas-Fermi 
radius in the case of a strongly interacting BEC. For the following discussion 
let L2 = L^{D, C) and H = H^'^D, C). 

For three dimensions, the Sobolev embedding states that 

Theorem 1. if^'^(D,R,) is continuously embedded in L^(D, R) if D is open 
and satisfies the cone condition and 2 < q < 6. 

Furthermore, one can obtain that H^'^{D, C) is continuously embedded in 
L'^{D,C) if D is open and satisfies the cone condition and q satisfies the 
inequality for the above theorem. For the one- and two-dimensional case this 
theorem also holds, but a much stronger result is available. See |25] for details 
on Sobolev spaces and the Sobolev embedding theorem. Next, we present our 
result for uniform convexity. 

Lemma 1. If we assume that the trapping potential is bounded away from 
zero, E{u) is uniformly and strictly convex. In other words, for all e > with 
Vtrapix) > e, there exists 5 > Q so that for all u and h E H 

E"iu)ih,h)>5\\h\\l. 



Proof. To see this we do a direct computation to obtain 

E"{u){h, h) = {\Vh{x)\^ + 2Vtrap{x)\h{x)\^) dx 

+ I 2g\h{x)\'^\u{x)\^ + Ag {^{h{x),u{x))fdx 

JD 

> l^{\Vhix)\'' + 2Vtrapix)\hix)\^) dx 

> 1^ {5\Vhix)\^ + 5\hix)\'') dx = 5\\h\\l 

where 5 = min(2e, 1). □ 

Furthermore, note that this assumption on the trapping potential does not 
lead to a loss of generality as given by the following lemma. 
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Lemma 2. Let E{u) he as in ([T]) and 



EJu) = E(u) + ef \u(x)\'^dx. 

JD 

Then for ^{u) = \u{x)\'^dx and he null{(3'{u)), E'{u)h = iff E'^{u)h = 0. 
Proof. Note that 

E'^{u)h = E'{u)h + e(3'{u)h. 
Thus if h E nuU{P'{u)) we have that P'{u)h = and 

E'^{u)h = E'{u)h + e(3'{u)h = E'{u)h. 
Thus we have that E[{u)h = iff = for h e null{(3' {u)) . □ 

Thus if the minimum of the trapping potential is zero, we can add a multiple 
of |M(a;)pda:; to the energy to obtain E^{u) as given above. Additionally, 
we see that if we have u so that E'^{u)h = for all h G null {13' {u)) , then 
E'{u)h = also for all h G null{[3' {u)) and the desired result is achieved. 

1.2 Minimization of the energy functional with the Sobolev gradient 

In this section we present the method of Sobolev gradients which can be 
applied to the problem of finding the minimizer of the Gross-Pitaevskii en- 
ergy functional ([1]). We define the Sobolev gradient and subsequently discuss 
how the normalization constraint ([2]) can be incorporated in the gradient. 
We pursue this approach since we propose a direct minimization for which 
a traditional Lagrange multiplier method is not suited. The motivation and 
background for this minimization closely follows 

Since ii^ is a continuously twice differentiable function from to IR C C, the 
Riesz representation theorem provides the result that for u E H, there is a 
member of H which we denote by VhE{u) so that 

E'{u)h = {h,VHE{u))H (3) 
for all h E H. We denote V hE{u) to be the Sobolev gradient of E at u. 
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We wish to minimize E using the method of steepest descent with this gra- 
dient. However, we first need to incorporate the normahzation constraint into 
our formulation. For u & H, let 



(3{u) = [ \u{x)\Mx 



Jd 



(4) 



which means that we want to minimize the energy functional ([T]) subject to the 
constraint (3{u) = N. We propose a direct minimization method using steepest 
descent with a gradient obtained with respect to a Sobolev inner product. In 
order to incorporate the normalization constraint into our formulation, we 
cannot use a traditional Lagrange multiplier method to add a multiple of the 
constraint to the energy as this would make the resulting functional unbounded 
from below. 

However, in drawing motivation from the method of Lagrange multipliers, we 
see that if such a minimum is achieved at u then the gradient of the energy 
functional E at u and the gradient of /? at m are parallel, and additionally 
(3{u) = N. Thus we seek u E H so that P{u) = N and if {h,'V hP{u))h = 
for h & H, then {h, V hE{u))h = also. Note that this translates into finding 
u so that P{u) = N and E'{u)h = for all h G null {(3' {u)). The rest of 
this section is devoted to a formulation of how such a u can be found in the 
scope of steepest descent. The general setup on how to include constraints in 
the gradient is explained in detail in chapter "Boundary and Supplementary 
Conditions" of [2B], but the essential ingredients are reviewed in the next 
paragraph. 

As previously stated, for each u & H one can find a member of if, denoted 
by V_h--E(m), so that equation ([3]) holds for all h & H. Note that for each 
u & H, the nuUspace of f3'{u) is a closed linear subspace of H. Thus for each 
u E H, there is a projection from H onto the nullspace of l3'{u). We denote 
this projection by Pu in this paper. Let Uq E H so that P{uq) = N and define 
z : [0, oo) ^ H so that 



Here z'{t) denotes the Frechet derivative as defined above and we will show 
that z is defined for alH > in section [T31 Essentially, we project the Sobolev 
gradient of E at z{t) onto the nullspace of P'{z{t)) since this leads to a constant 
f3{z). This can be seen by considering 



where we used that Pz(t) is the projection of H onto the nullspace of P'{z(t)). 
Hence, f3{z) is constant and if m = lim^^oo then (3{u) = (5{uq). 



z(0) = Mo and z'(t) 



P,(^t)^HE{z{t)). 



(5) 



{m)\t) = /3'izit))z'it) = -f3'izmp<t)VHEizm = o 
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Thus by projecting the Sobolev gradient of E at z(t) into the nuUspace of 
f3'{z{t)) for each t, we obtain that z{t) satisfies the constraint for all t and in 
the limit u. We will also show that if m = linif^oo z{t) for z as given in equation 
(jS]), then E'{u)h = for all h in the nullspace of f3'{u) which will give us the 
u that we desired above. We also provide a convergence result for z in section 
11.41 In particular, if we assume that the trapping potential is bounded away 
from zero and non-negative we obtain that lim^^oo z{t) exists. 

In our formulation we will need to know what the projection is. The next 
part of this work is devoted to finding an expression that we can use both in 
our analysis and in doing our simulations. 



1.3 Properties of the projection onto the nuUspace of f3'{u) 

In this section we discuss the previously mentioned projection P„ which con- 
stitutes an essential part of the Sobolev gradient. Before we show the explicit 
representation for the projection P„, we need to introduce the following defi- 
nition and a theorem which is proved in [27] . 

Definition 2. Let f E L'^(D), and define af from H^''^{D) to R so that 
af{u) = (m, /)l2(£,). Since H^'^{D) is continuously embedded in L'^{D), af 
is continuous from H^''^{D) to R.. Thus there exists v G H^''^{D) so that 
{u,f)L^{D) = {u,v)h^,2(^d) for allu G H^''^{D). Define Mf = v. 

Theorem 2. M is mjective. M G L{X,Y) with X = H^^'^{D) or L'^{D) and 
Y = H^'^{D) or L?'{D). In each case, the operator norm of M is less than or 
equal to one. 

These two results yield the following relationship between the two inner prod- 
ucts (-, ■)l2(£,) and (-, ■)m,2^D) 

{uJ)lHD) = {u,Mf)Hi.2(^D) 

for all u G H^'^{D) and / G L'^{D). 

We also have the following theorem for which we start with / G H^'^{D) and 
let W f = V/. Thus, W is a closed densely defined linear operator from L'^{D) 
to [L'^{D)Y whose adjoint we denote as W*. The following representation for 
M is proved in [28] . 

Theorem 3. Let M be the transformation defined in definition\^ then M = 

{i + w*w)-\ 
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Using the same formulation given in Definition [2] for L'^{D, C) and H^''^{D, C), 
one can also obtain an operator : L'^{D, C) — >• H^''^{D, C) so that {h, u) ^2 = 
(/i, M(sju)h for all /i G if and u & L'^. Then one has the following result. 

Theorem 4. If u = r + is e L'^{D, C), t/ien Mcm = Mr + iMs. 

Proof. We have that {h,u)L2 = {h,M^u)H for all h E H. First suppose that 
h = f + Oi which yields 

^{h, u)l2 = ^{h, Mcu)h = if, ^{Mcu))hi,^d) , 

where 3? denotes the real part. Furthermore, we observe that 

^{{h,u)L2) = {f,r)L2(D) = {f,Mr)Hh2(D)- 

and thus 

{f,Mr)H^.^D) = {fMMcu))m.2^D) 
for all / G H^'^{D) and it must be that Mr = 3?(Mc-u). A similar argument 
shows that Ms = 53(Mc'u), where ^3 denotes the imaginary part. □ 

From this result it follows that has the same properties as M. From now on 
we will not distinguish between M and in our notation. For u G L'^{D, C), 
Mu = Mr + iMs. 

Using these definitions and results, we obtain an explicit formula for the pro- 
jection Pu 



Puh = {I- Ql{QuQl)-'Qu)h = h- J^^'l^^f Mu , (6) 

M.{U, Mu) 1,2 

which is derived in Appendix \M 

Looking at this expression for the projection P^, we note that as u varies, 
the associated projections vary in a continuous way. We will use this result 
to obtain global existence for z. In more precise terms, we have the following 
result which is proved in Appendix [Bl 

Proposition 1. Suppose {«„}„>! is a sequence of members of H that con- 
verges in L"^ to u ^ & H . Then the sequence {Pu„}n>i converges in L{H, H) 
to Pu- Furthermore there is a constant m so that ||-Pu„ — -P«|| < "Till-u^ — u\\i2. 
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1.4 Global existence, uniqueness and convergence of the Sobolev gradient 



In this section we first show that z as given in equation ^ is defined for all 
t > by using Proposition [TJ Additionally, we prove that linii^oo z{t) exists 
in H if the trapping potential is bounded away from zero by some positive 
number e. 

Theorem 5. z as given in ([5]) is uniquely defined for allt>0. 

Proof. Since E is continuously twice differentiable, the map u — *• VhE{u) 
from H to H is a Lipschitz map. Furthermore, due to the result obtained in 
Proposition [H we have that the map u —>■ PyV hE{u) is also a Lipschitz map. 
This implies that we have local existence for equation Suppose that there 
exists a number T so that z as defined in equation ([5]) can only be defined for 
< t < T. To see that we have global existence, let < a < 6 < T and note 
that 

\\z{b)-z{a)\\l<(^j\\z'{t)\\ndt 

^ \\P,(^VHE{z{t))\\Hm 

<{b-a) ['\\P,^t)VHE{z{tmj,dt 
<T f'\\P,(^t)VHE{zm\ldt. 

J a 

Furthermore, we have 

{E{z))'{t) = E'{z{t))z'{t) 

= {z'{t),VHE{z{t)))H 

= -{P4t)VHE{z{t)), VHE{z{t)))H 

= -{P\,)VHE{z{t)),VHE{z{t)))H 

= -{P.it)VHE{z{t)),P,f^)VHE{z{t)))H 

= -\\P,^t)yHE{z{t))\?H<^ 

for all t < T which implies that E{z) is decreasing on [0,T). Additionally, 

So \\Pz{t)^HE{z{t))\\ldt < oo since 

/' \\P,(t)VHE{z{t))fH dt = E{z{a)) - E{z{b)) < E{z{0)) . 

J a 
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Thus limf^r- z{t) = z{T) exists and one can therefore define z on [0,T + e) 
for some positive number e, which means that z{t) must exist for all t > 0. 
Uniqueness follows from basic existence and uniqueness for ordinary differen- 
tial equations. □ 



This global existence and uniqueness result for z as given in equation (jS]) is 
accompanied by global convergence to an element of the Sobolev space u E H. 
Our result is an adaptation of what is presented in [26] and implies strong 
convergence in the H norm. This is a stronger result than that of using a 
minimizing sequence to obtain weak H convergence or strong convergence 
as has been previously done. For our proof, suppose that the trapping potential 
Vfrap is bounded away from zero by some positive number e. This does not lead 
to any loss of generality as previously described. Recall that this assumption 
implies that there is a positive number 6 so that 

E"{u){h,h) > 6\\h\\]j. 

Theorem 6. Suppose that z is given by equation dSD, with V hE{uq) ^ 0, 
then 

lim z{t) = u 

t—KX> 

exists. Furthermore, there exist constants m andc so that \\u — z{t)\\H < me~^*, 
and E'{u)h = for all h E null{l3'{u)). 

Proof. Let g{t) = E{z{t)) then we know already that 

g'it) = -\\P,^,)WHEizit))fH- 

Moreover, let G : ^ if be defined by G{x) = P^S/ hE{x) then G is and 
thus 



g"{t) = -2{G\z{t))z\t),G{z{t)))H . 
Now note that ii h E null{j3'{z{t)), then 



E'izit))ih) = {h, WHE{z{t)))H = {P.it)h, WHE{z{t)))H = {h, ^(^(t)))^ 
and therefore 



E"izmKz'it)) = {h,G'izit))z'it))H . 
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Since 

z'{t) = -P,^t)VHE{z{t)) = -G{z{t)) 
is in the nullspace of P'{z{t)), one has 

E"{zmAt),^\t)) = {z'{t),GUt))z'{t))H 

= {-G{z{t)\G'{z{t))z\t))n = ^-^ 

Making use of Lemma [1] we consequently have 



and hence 
9"{t) 



> 26. 



g'it) 

Integrating both sides from to t, we get that 

-\n{-g'{t)) + \n{-g'm>25t 
and thus there is a constant m so that 

< -g'it) < for all t. 

More specifically m = \\Pz(o)V hE{z{0))\\1[. From this we see that 



71+1 \ 2 rn+l rn+1 rn-\-l 

Ik'H^dt </ ||/||^dt = -/ g'dt<m e-^^'dt 
and therefore 



oo 

\\z'\\h dt < go 



which implies that lim^^oo z{t) = u exists. The rate of convergence is given by 



■Sn 



\u-z{n)\\H< I - )e~ 
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To see that E'{u)h = for all h G null{(3'{u)), recall that by Proposition [H 
the map u ^ from H to L{H, H) is Lipschitz. 



Furthermore since E is C^, the map u —* VhE{u) is Lipschitz also from H to 
H. Thus the map u — > PyV hE{u) is Lipschitz and from this it follows that 

PuVhE{u) = lim P,^t)VHE{z{t)). 

t — *oo 

Since 

roo rco 

/ \\z'\\Hdt= / \\P,^t)VHE{z{t))\\Hdt 
Jo Jo 

is finite, then it must be that PyV hE{u) = 0. Now let h G null{(3' {u)) then 

E'{u)h = {h, VhE{u))h = {Puh, VhE{u))h = {h, PuVhE{u))h = , 
which concludes the proof. 

□ 
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2 Numerical results 



In the second part of our contribution we want to apply tlie results that 
we obtained in the detailed discussion about the theoretical background of 
the Sobolev method and its application to calculate the ground state of the 
Gross-Pitaevskii energy functional. In particular we want to demonstrate the 
extraordinary convergence properties of this method and highlight its numer- 
ical efficiency. For this purpose we do not use the commonly known harmonic 
oscillator, but the more challenging generalized Mexican hat potential as our 
trapping potential which in natural units reads 



for the (i-dimensional case. For our numerical simulations we choose A = j^, 
B = 16 and - to break the spherical symmetry - Ci = 1.0, C2 = 1.5 and 



Furthermore, we vary the coupling constant g from 10~^ over 10° up to 10^ 
such that we cover the regime from a nearly interactionless gas of bosons to 
one which is interaction dominated. For g = 10~^ the bosons are only weakly 
interacting, whereas they are strongly interacting for g = 10^. For a numerical 
realization of the minimization we also need the particle number N which 
determines the normalization constraint. The particle numbers we choose for 
our illustration purposes are = 10^ for the one- dimensional case, = 10^ 
in two dimensions and = 10^ in the three-dimensional case. 

Before we discuss our results, a few words are in place about our numerical 
implementation of the method of Sobolev gradients. The details about our 
implementation are summarized in Appendix C and we briefly want to mention 
that we use a discretization in position space and an implementation of the 
differentiation by means of central differencing [29]. Furthermore, we use an 
Euler iteration to discretize equation ([5]) and optimize the step size for each 
iteration. Once the relative change in energy from one iteration step to the 
next is small enough - in our case this means less than 10~^ for ID, 10~^ 
for 2D and 10""^ for 3D due to the speciflc choice of the potential and the 
respective particle numbers - we switch to the well-known Newton method in 
order to flnd the exact solution of the stationary Gross-Pitaevskii equation 

~V^u{x) + Vtrap{x)u{x) + g\u{x)\^u{x) = fiu{x) , 
where n denotes the so-called chemical potential. 




C3 = 2.0. 
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The ground state of the Gross-Pitaevskii energy functional has to fulfill this 
equation and only few steps with the Newton method are necessary to arrive 
at a desired accuracy of 10"^ in each component. 

To clearly demonstrate the superior convergence properties of our numerical 
simulations with the Sobolev gradient, we choose the initial wave function to 
be a real- valued random field with zero boundary conditions which on average 
is as far as possible from the desired solution. All simulations were run on a 
Dell PowerEdge 2900 Server with four Intel Xeon 5160 CPUs at a frequency 
of 3.00 GHz and a total of 16 GB of RAM. The MATLAB program for our 
simulations is a 64-bit Linux version 7.5.0.338 (R2007b). 



2.1 Simulations in ID 



AT, 


9 




#S,Min 


#S,Max 


#N,Min 


#N,Max 




tMax 




10-1 


4.8938 


130 


241 


17 


18 


5.2-10-1 


8.6-10-1 


27 


10° 


19.831 


98 


123 


21 


23 


4.4-10-1 


5.9-10-1 




lOi 


92.607 


54 


86 


19 


23 


4.2-10-1 


5.7-10-1 




10-1 


4.8952 


156 


278 


18 


19 


6.6-10-1 


1.1-10° 


28 


10° 


19.831 


109 


179 


20 


23 


5.0-10-1 


7.4-10-1 




lOi 


92.607 


80 


121 


20 


20 


4.1-10-1 


5.5-10-1 




10-1 


4.8956 


153 


309 


18 


20 


7.6-10-1 


1.4-10° 


29 


10° 


19.832 


124 


208 


20 


23 


6.5-10-1 


1.0-10° 




lOi 


92.608 


93 


147 


21 


21 


5.3-10-1 


7.6-10-1 



Table 1 

Simulation results with the Sobolev gradient for the ground state of the ID Gross- 
Pitaevskii energy functional with a Mexican hat potential. 



In table [2TT] we see the results of the one-dimensional simulations for a particle 
number of = 10^ and a grid length of = 10. For various discretization 
numbers and coupling strengths g we run our simulations to obtain the 
ground state for 10 different random initial fields. We show the respective 
chemical potential fi in dimensionless units, the minimal and maximal num- 
ber of iteration steps with the Sobolev gradient #s,Min/Max, the minimal and 
maximal number of iteration steps with the Newton method #N,Min/Max and 
the minimal and maximal CPU time tMin/Max for the simulation in seconds. 
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Fig. 1. Density distribution = |n(x)p of the ground state of the ID Gross-Pi- 
taevskii energy functional with a Mexican hat potential for a particle number of 
N = 10^. The number of grid points for the simulation is Nx = 2^ and the grid size 
is Lx = 10. The interaction strength is g = 10~^ in subplot (a), g = 10*^ in subplot 
(b) and g = 10^ in subplot (c). 

which was obtained in our MATLAB simulation with the built-in functions 
tic and toe. 

We can clearly observe the dependence of the number of iteration steps with 
the Sobolev gradient on the different random initial fields. In general the 
number of necessary steps decreases for an increasing coupling constant. In 
contrast, the number of iterations with the Newton method is almost inde- 
pendent of the random initial field and shows a negligible dependence on the 
coupling constant. Furthermore, it is worth mentioning that the simulation 
time is independent of the number of grid points since the complexity of the 
minimization in one dimension is fairly low and thus the computational over- 
head is dominant. 

The shape of the density distribution Ux = \u{x)\'^ is depicted in figure [H for 
Nx = 2^ and interaction strengths of g = 10~^, 10° and 10^. We can notice the 
change from a clear double Gaussian shape for a small interaction strength 
towards a strongly interaction broadened shape at a large interaction strength. 

2.2 Simulations in 2D 

In table 12.21 we compare the parameters for the two-dimensional simulation 
for a particle number of = 10^ and a grid length of Lx = Ly = 10. The 
discretization for both spatial parameters x and y was chosen to be equal by 
setting Nx = Ny. As previously, we show the respective chemical potential fi in 
dimensionless units, the minimal and maximal number of iteration steps with 
the Sobolev gradient #s,Min/Max; the minimal and maximal number of iteration 
steps with the Newton method #N,Min/Max and the minimal and maximal CPU 
time tMin/Max for the simulation in seconds. 

For the two-dimensional case the number of iteration steps with the Sobolev 
gradient also depends on the different random initial fields as was the case 
for the one- dimensional simulations. The number of iteration steps with the 
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9 




#S,Min 


#S,Max 


#N,Min 


#N,Max 




^Max 




10-1 


5.6656 


350 


374 


21 


23 


5.7-10° 


9.9 - 10° 


26 


10° 


23.557 


179 


190 


23 


26 


5.1 - 10° 


5.4 - 10° 




lOi 


124.44 


87 


103 


25 


27 


4.4 - 10° 


5.3 - 10° 




10-1 


5.6836 


427 


451 


20 


21 


1.1 - 10^ 


1.5 - 102 


2^ 


10° 


23.564 


152 


243 


23 


23 


4.6- lOi 


6.1 - IQi 




lOi 


124.44 


93 


176 


25 


30 


3.5 - lOi 


6.3 - IQi 




10-1 


5.6883 


615 


990 


24 


24 


9.2 - 10^ 


2.3 - 10^ 


28 


10° 


23.566 


324 


626 


22 


24 


6.0- 10^ 


1.1 - 10^ 




lOi 


124.44 


142 


236 


26 


28 


4.6- 10^ 


6.7- 10^ 



Table 2 

Simulation results with the Sobolev gradient for the ground state of the 2D Gross- 
Pitaevskii energy functional with a Mexican hat potential. 



Newton method has again a negligible dependence on the different random ini- 
tial fields. Overall, the number of necessary steps has a strong dependence on 
the coupling constant and is significantly lower for a large coupling constant, 
which is due to a much smoother behavior of the respective ground state wave 
function. Now the simulation time depends on the number of grid points since 
the complexity of the minimization in two dimensions is rapidly growing. Com- 
bining the increasing number of iteration steps for twice as many grid points, 
the simulation time approximately grows by a factor of four which is exactly 
what one would expect since the complexity of a two dimensional system also 
grows by this factor when the number of grid points is doubled. Nevertheless, 
the simulation time is on the order of several minutes despite a maximum 
system size of 2^^ grid points. This clearly demonstrates the efficiency of the 
method of Sobolev gradients. 

To illustrate the shape of the density distribution rix = \ax\'^, we depict in 
figure [2] cuts along the x-axis for y = (black lines) and along the y-axis 
for a; = (blue lines). For these plots we used a discretization number of 
= 2^ and interaction strengths of g = 10~^, 10° and IQi. Once more we can 
notice the change from a clear double Gaussian shape for a small interaction 
strength towards a strongly interaction broadened shape at a large interaction 
strength. Since we deal with an anisotropic trapping potential the width in 
the ^/-direction is smaller than in the x-direction. 
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Fig. 2. Cuts through the density distribution Ux = |u(a;)p of the ground state of the 
2D Gross-Pitaevskii energy functional with a Mexican hat potential. The cuts are 
along the x-axis for y = (black lines) and along the y-axis for x = (blue lines). 
The results are for a particle number of = 10'^, the number of grid points for the 



simulation is = 
strength is g = 10" 
(c). 



Ny = 2 and the grid size is Lx = Ly = 10. The interaction 
in subplot (a), 5 = 10° in subplot (b) and g = 10^ in subplot 



1 ; 



2.3 Simulations in 3D 



Nx 


9 




#S,Min 


#S,Max 


#N,Min 


#N,Max 




^Max 




10-1 


9.8716 


767 


955 


24 


28 


1.2 • 10^ 


2.5 • 10^ 


25 


10° 


45.302 


297 


341 


23 


28 


5.1 • 10^ 


8.2 • 10^ 




lOi 


218.23 


133 


152 


25 


29 


2.5 • 10^ 


4.8 • 10^ 




10-1 


9.9784 


755 


987 


24 


29 


1.7- 10"^ 


2.4 • 10^ 


26 


10° 


45.414 


354 


421 


22 


23 


7.8 • 10^ 


1.1 • 10^ 




lOi 


218.18 


199 


215 


25 


29 


2.9- 10^ 


5.7- 10^ 




10-1 


10.011 


1078 


1485 


23 


25 


1.5 • 10^ 


2.3 • 10^ 


2^ 


10° 


45.443 


640 


769 


23 


23 


7.8 • 10"^ 


1.4 • 10^ 




lOi 


218.21 


371 


466 


26 


33 


5.0 • lO'' 


6.9 • 10^ 



Table 3 

Simulation results with the Sobolev gradient for the ground state of the 3D Gross- 
Pitaevskii energy functional with a Mexican hat potential. 



Eventually we compare the parameters for the three-dimensional simulation 
in table 12.31 We used a particle number of = 10^ and a grid length of 
Lx = Ly = Lz = 10. The discretization for all spatial parameters is the same 
and thus A^^ = Ny = N^. Again we show the respective chemical potential n in 
dimensionless units, the minimal and maximal number of iteration steps with 
the Sobolev gradient #s,Min/Max and with the Newton method #N,Min/Max and 
the minimal and maximal CPU time tMin/Max for the simulation in seconds. 
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Fig. 3. Cuts through the density distribution Ux = \u{x)\'^ of the ground state of 
the 3D Gross-Pitaevskii equation with a Mexican hat potential. The cuts are along 
the X-axis for y,z = (black lines), along the y-axis for x,z = (blue lines) and 
along the z-axis for x,y = (red lines). The results are for a particle number of 
N = 10^, the number of grid points for the simulation is Nx = Ny = = 2^ and 
the grid size is Lx = Ly = Lz = 10. The interaction strength is g = 10~^ in subplot 
(a), 5f = 10'^ in subplot (b) and g = 10^ in subplot (c). 

As before the number of iteration steps with the Sobolev gradient depends on 
the different random initial fields whereas they do not influence the number of 
iteration steps with the Newton method. The number of necessary steps also 
strongly depends on the coupling constant and is lower for a large coupling 
constant, where the respective ground state wave function has a smoother be- 
havior. The simulation time is highly dependent on the number of grid points 
due to the complexity of the three dimensional simulations. The simulation 
time is now on the order of one day for a maximum system size of 2^^ grid 
points. However we want to point out, that this is not a standard computation 
time for a given potential but represents the upper limit for a ground state 
simulation with our method, since we started as far as possible from the final 
solution and with the most unfavorable initial state. 

As before, the form of the density distribution Hx = \u{x)\'^ again changes 
from a clear double Gaussian shape for a small interaction strength towards 
a strongly interaction broadened shape at a large interaction strength. Since 
we deal with an anisotropic trapping potential the width in the z-direction is 
smaller than in the ?/-direction, which in turn is dominated by the x-direction. 
In figure [3] we illustrate this behavior by depicting cuts through the density 
distribution along the x-axis for y,z = (black lines), along the y-axis for 
x,z = (blue lines) and along the z-axis for x,y = (red lines). The plots 
are for Nx = 2'^ and interaction strengths oi g = 10^^, lO'' and 10^. 
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Conclusion 



In this paper we performed a detailed study of the minimization of the Gross- 
Pitaevskii functional using the method of Sobolev gradients and the trajectory 
given in equation ([5]). In the analytical part of our work, we obtained global 
existence and uniqueness for this trajectory as well as global convergence to 
a minimizer of the Gross-Pitaevskii functional. Furthermore, in our numerical 
part we used finite differences to discretize this trajectory and were able to 
find stationary solutions in one, two, and three dimensions. In our study we 
found that the advantages our method presents are high numerical stability, 
fast convergence to the desired solution, and versatility in handling various 
parameters such as the trapping potential, the coupling constant and initial 
estimates. The main new contribution of the analysis we performed is that we 
give an explicit minimizing sequence that converges to the minimizer of the 
Gross-Pitaevskii functional. In contrast to previous comparisons of numeri- 
cal methods for the minimization of this functional, our goal is not to find 
the smallest computation time for a given potential with well known starting 
point. Here we show that using an initial random field, that is without know- 
ing anything about the final ground state, we arrive at the desired solution 
of the stationary Gross-Pitaevskii equation without any additional technical 
tricks. This is due to the global convergence of our minimizing sequence and of 
particular interest for arbitrary external potentials that are not as well known 
as the standard choices, like a harmonic potential. 

The work we present in this paper opens the door for many interesting studies 
of the Gross-Pitaevskii energy and equation. As a follow up numerical project, 
we are currently studying the performance of this method for the case of the 
Gross-Pitaevskii energy with rotation, which corresponds to a BEG in a rotat- 
ing frame. We plan on numerically investigating the formation of vortices, but 
will also consider analytical studies about existence and convergence proper- 
ties of vortex lattices. On the other hand, we are currently working on the 
adaption of our scheme to study the time dependent Gross-Pitaevskii equa- 
tion. General ideas of how such an investigation can be done are presented in 

m- 

Apart from physical applications our method also provides a starting point 
for mathematical investigations, in particular from the perspective of nonlin- 
ear semigroup theory. Defining the operator T by T{t)x = z{t) where z is 
given in with z{0) = x, we note that T is a strongly continuous nonlinear 
semigroup. We propose to look for the Lie generator of this group as well 
as studying properties of it. Since the Sobolev gradient provides an explicit 
construction of the trajectory with global convergence, we expect various in- 
teresting properties of this generator. 
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A Explicit representation of the projection P„ 

In this section we first provide the general setting which is used to derive a 
formula for P„ and thereafter derive the explicit representation of P^. A general 
result from functional analysis states that if Q is a linear transformation from 
one Hilbert space X to another Hilbert space Y and if Q*, the adjoint of Q, 
as a continuous linear transformation from F to X has closed range, then X 
has a unique decomposition as A = null{Q) © range{Q*). Furthermore if P 
is the orthogonal projection of X onto null{Q), then / — P is the orthogonal 
projection of X onto the range{Q*) . Observe that Q*{QQ*)~^Q is symmetric 
from X to X, idempotent, has range in range{Q*), and is fixed on this set. 
Thus Q*{QQ*)~^Q is the orthogonal projection of X onto range{Q*). This 
makes I — Q*{QQ*)~^Q the orthogonal projection of X onto null{Q). 

Now we apply this general setting to our case with Qu = P'{u)- Then for 
heH = H^'\D, C) we obtain 

Quh = (3'{u)h = 2^{h,u)L2. 

We now want to compute the adjoint of Qu as a continuous linear transfor- 
mation from i7 to R C C. Recall the definition of M from Definition [21 For 
r G R and h e H 

{Quh, r)R = Quh *r = 2^{h, u)L2*r = {h, 2rMu)H- 

Thus we see that Q*r = 2rMu. Now suppose the sequence ?/„ = Q*r„ = 
2rnMu is in the range of and converges to f G if. Then this sequence is 
Cauchy and for e > 0, there is N so that if m,n > N 

\\2{rn-rra)Mu\\H = 2\rn - rm\\\Mu\\H < e. 
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This implies that {r„}n>i is a Cauchy sequence in R and hence converges to 
some r G R. Thus {2r„M-u}„>i converges in H to 2rMu and v = 2rMu is 
in the range of Q*, which imphes that the range of is closed. Thus the 
orthogonal projection of H onto the nuUspace of Qu is given by 



Pu = i-q:{Quq:)-'Qu. 

Now observe that QuQuf — Qui'^fMu) and since Quh — 23?(/i, u)l2 for h e H, 
one has 



Qu{2rMu) = 2^{u, 2rMu)L^ = m{u, Mu)l^ * r . 

Thus QuQIt = 4:'Si{u, Mu) L2 * r and {QuQD^^ exists if m 7^ 0. Furthermore, 
{QuQu)~^ is continuous and given by 



Consequently we arrive at 



Q:{QuQ:)-'Quh^Q:{QuQ:r'2^{u,h)L^ 

and obtain the explicit representation of the projection P„ as 

Puh = (/ - QliQuQir'Qu) h = h- ^^§£-Mu . 



B Lipschitz property of the projection P„ 



Proposition 1. Suppose {un}n>i is a sequence of m,emhers of H that con- 
verges in K to u ^ ^ H . Then the sequence {P„„}n>i converges in L{H, H) 
to Pu- Furthermore there is a constant m so that \\Pun ~ -Pull ^i^\\un — u\\l2. 

Proof. The proof will be given in two steps. Let h E H with \\h\\H = 1- First 
note that 



{Mu,u)l2 = {Mu,Mu)h = \\Mu\\l, 
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which yields that Mu 7^ since « 7^ and M is injective. Furthermore, we 
recall that M E L{X,Y) where X = H,L'^ and Y = H,L'^. To minimize 
notation we write that ^{f,g)L^ = {f,g)K- 

Thus we have 



{Un, h)K _ {u,h)K _ {Un, h) K\\Mu\\jj - {u, h) K\\MUn\\jj 



MUnWn \\Mu\\jj \\Mu\\jj\\MUn\\H 

~ II M7/JI?rll M7/,„ll?r 



< U ^ 



||MM|||^||MM„||f^ 

II Mm 11^ II Mm„ 111^ 



^ 11^ - ^n ||L2 



< 11^ - '"nllLHI'^llia + II^IU2|(||Mm||^ - ||. 

~ mi 
for some number mi. Additionally, we note that 



\\MU\\j^\\MUn\\H 

\Wl^ + ||^IU ^|(||Mm||^-||Mm.||^)| 



||Mm||^- ||Mm„||^ ^\{Mu,u)k- {Mu^,u.^)k\ 



< \{Mu,U- Un)K\ + \{MU - MUn,Un)K\ 

< ||Mm||l2||m - m„||l2 + ||M('u„ - 'u)||L2||ti„||i2 

<m2\\Un - u\\l2 



for some number 1112 ■ Now if we let 



{Un, MUn)K 



and c 



{'l.h)K 

{u, Mu)k ' 



it is clear that there is a constant ma so that \cn — c\ < msjlwn — m||l2- Thus 
we obtain 



\\{Pun - Pu)h\\H=\\CnMUn - cMu\\h 

< ^CnMUn — CnMu^H + jjCnMli — cMu^h 

< |C„|||MK - u)\\h + \Cn - c\\\Mu\\h 

and consequently there is a constant m so that 



||(-P«„ - Pu)h\\H<m\\u-Un\\L'^ , 
which concludes the proof. 



□ 
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C Numerical implementation 



We discretize the continuous problem in position space and for convenience 
we only consider the discretization of a one-dimensional problem, since the 
extension to higher dimensions is self-explanatory. The discrete form of the 
spatial variable x, restricted to the finite interval [—L^, L^], reads 



X = (xi, ...,xnJ where x„ = 2L^ - - j , ^ < n < N^, e¥i . 

The length has to be sufficiently large such that the wave function that 
we are interested in is numerically zero outside the interval [—L^jL^]. The 
existence of such a length is guaranteed by our requirement that the exter- 
nal potential diverges for \x\ — > oo. The number of grid points is a very 
crucial parameter for the simulation because it is related to the grid spacing 
= 2Lx/Nx and thus determines the possible resolution of the scalar func- 
tion u. The wave function u on the grid introduced above is also represented 
as a vector 



u = (til, . . . , un^) where Un = u{xn), I < n < 



We approximate the spatial derivatives of the wave function via first order cen- 
tral differencing, which provides an accuracy for the first derivative approxima- 
tion at cell centers of order Aj^. Denoting the cell centers by Cj = (a^i+i +Xi)/2 
we obtain 



fi+l ~ fi 



For this scheme the discretized version of the operator W from the analytical 
part of this contribution is an A^^ — 1 x A^^ matrix that reads 



W ^ — 



1 

A. 



-1 1 

0-1 1 ... 

... 0-1 10 
0-11 



such that W{f) 



( h-h \ 



Ax 



It is important to recall that this approximation converges towards the exact 
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derivative in the limit of an infinitely fine grid. We use different values for the 
number of grid points A^^; which allows us to show that the Sobolev gradient 
indeed converges. 

Using this differencing, we sec that the discrete inner product for H^'^{D), 
denoted by (•, •)5 is given by the following. For /, g being valued functions 

{f:9)s^{f:g)N+{W{f),W{g))^^, 

where {■,-)n denotes the inner product. Note that 

(/, 9)s = if: 9)n + (/, W*W{g))i, = (/, (/ + WW) {g))^. 

I + Vr*Vr is positive definite and hence injective. Furthermore, it is invertible 
which allows us to obtain a result analogous to the infinite dimensional case 
for the relationship between the H^'^ inner product and the inner product. 
We see that 

(/, {I + W*W)-'g)s = {f.g)N 
and thus the analogous finite dimensional M as given in the first part is 

{i + w*w)-\ 

Now, all we need for the numerical simulation is an easily accessible procedure 
to calculate the Sobolev gradient. We restrict ourselves to real- valued functions 
M, since it can be shown, that the ground state can be chosen to be real-valued. 
The Frechet derivative of the energy functional then reads 

E'{u)h = {^lyu{x), Vh(x)) + 2Vtrap(x)n{u(x), h{x)) 
+ 2g\u{x)\^^{u{x), h{x))) dx . 

where (/, g) is a short-hand notation for the product of / and the complex 
conjugate of g. By adding and subtracting 'R{u{x), h{x)) da; from this term, 
one obtains 

E'{u)h = 3?(/i, u)h + 2VtrapU + 2g\u\^u - u)l2 . 

Making use of the previously mentioned relationship between the if^'^ inner 
product and the LP' inner product we arrive at 

E'{u)h = 3fJ(/i, u + M{2VtrapU + 2g\u\\ - u))h 
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where M is as in Definiton [2l Recall that the Sobolev gradient of E at m, 
VhE{u), was defined to be the element of H so that 

E'{u)h = {h,VHE{u))H = ^{h,VHEiu))H 
as E'{u)h is real valued. Therefore, the Sobolev gradient of ii^ at u is given by 

VhE{u) =u + M{2VtrapU + 2g\u\'^u - u) . 
Incorporating the projection of H onto the nullspace of I3'{u) one sees that 



As soon as the minimization of the energy functional has converged in the 
sense of the relative change in energy from one iteration step to the next is 
less than 10"'^, we switch to the well-known Newton method in order to find 
the exact solution of the stationary Gross-Pitaevskii equation 



with the chemical potential potential jj,. Hence it only remains to show, how 
the calculation of the chemical potential is implemented. For this purpose, 
we multiply both sides of this equation by the complex conjugate of u and 
integrate over the domain D. Thus, we arrive at 



PuVhE{u) = VhE{u) 



^{u,VhE{u))l2 



Mu. 



V'^u{x) + Vtrap{x)u{x) + g\u{x)\^u{x) = fiu{x) 
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